Water Surface Elevation in 3D

import xarray as xr

from schimpy import schism_mesh

import pandas as pd
import holoviews as hv
from holoviews import opts, dim

Open output dataset containing the water surface elevations

ds = xr.open_mfdataset('../tests/data/m1_hello_schism/outputs/out2d_*.nc', concat_dim='time', combine="nested",
                  data_vars='minimal', coords='minimal', compat='override')
ds

<xarray.Dataset>
Dimensions:                  (time: 144, one: 1, nSCHISM_hgrid_node: 2639,
                              nSCHISM_hgrid_face: 4636,
                              nSCHISM_hgrid_edge: 7274,
                              nMaxSCHISM_hgrid_face_nodes: 4, two: 2)
Coordinates:
  * time                     (time) datetime64[ns] 1999-12-31T16:20:00 ... 20...
    SCHISM_hgrid_node_x      (nSCHISM_hgrid_node) float64 dask.array<chunksize=(2639,), meta=np.ndarray>
    SCHISM_hgrid_node_y      (nSCHISM_hgrid_node) float64 dask.array<chunksize=(2639,), meta=np.ndarray>
    SCHISM_hgrid_face_x      (nSCHISM_hgrid_face) float64 dask.array<chunksize=(4636,), meta=np.ndarray>
    SCHISM_hgrid_face_y      (nSCHISM_hgrid_face) float64 dask.array<chunksize=(4636,), meta=np.ndarray>
    SCHISM_hgrid_edge_x      (nSCHISM_hgrid_edge) float64 dask.array<chunksize=(7274,), meta=np.ndarray>
    SCHISM_hgrid_edge_y      (nSCHISM_hgrid_edge) float64 dask.array<chunksize=(7274,), meta=np.ndarray>
Dimensions without coordinates: one, nSCHISM_hgrid_node, nSCHISM_hgrid_face,
                                nSCHISM_hgrid_edge,
                                nMaxSCHISM_hgrid_face_nodes, two
Data variables:
    minimum_depth            (one) float64 dask.array<chunksize=(1,), meta=np.ndarray>
    SCHISM_hgrid             (one) |S1 dask.array<chunksize=(1,), meta=np.ndarray>
    depth                    (nSCHISM_hgrid_node) float32 dask.array<chunksize=(2639,), meta=np.ndarray>
    bottom_index_node        (nSCHISM_hgrid_node) int32 dask.array<chunksize=(2639,), meta=np.ndarray>
    SCHISM_hgrid_face_nodes  (nSCHISM_hgrid_face, nMaxSCHISM_hgrid_face_nodes) float64 dask.array<chunksize=(4636, 4), meta=np.ndarray>
    SCHISM_hgrid_edge_nodes  (nSCHISM_hgrid_edge, two) float64 dask.array<chunksize=(7274, 2), meta=np.ndarray>
    dryFlagNode              (time, nSCHISM_hgrid_node) float32 dask.array<chunksize=(72, 2639), meta=np.ndarray>
    elevation                (time, nSCHISM_hgrid_node) float32 dask.array<chunksize=(72, 2639), meta=np.ndarray>
    depthAverageVelX         (time, nSCHISM_hgrid_node) float32 dask.array<chunksize=(72, 2639), meta=np.ndarray>
    depthAverageVelY         (time, nSCHISM_hgrid_node) float32 dask.array<chunksize=(72, 2639), meta=np.ndarray>
    dryFlagElement           (time, nSCHISM_hgrid_face) float32 dask.array<chunksize=(72, 4636), meta=np.ndarray>
    dryFlagSide              (time, nSCHISM_hgrid_edge) float32 dask.array<chunksize=(72, 7274), meta=np.ndarray>

xarray.Dataset

• Dimensions:
  • time: 144
  • one: 1
  • nSCHISM_hgrid_node: 2639
  • nSCHISM_hgrid_face: 4636
  • nSCHISM_hgrid_edge: 7274
  • nMaxSCHISM_hgrid_face_nodes: 4
  • two: 2
• Coordinates: (7)
  • time
    (time)
    datetime64[ns]
    1999-12-31T16:20:00 ... 2000-01-...

    i23d :

    0

    base_date :

    2000 1 1 0.00 8.00

    standard_name :

    time

    axis :

    T

    array(['1999-12-31T16:20:00.000000000', '1999-12-31T16:40:00.000000000',
           '1999-12-31T17:00:00.000000000', '1999-12-31T17:20:00.000000000',
           '1999-12-31T17:40:00.000000000', '1999-12-31T18:00:00.000000000',
           '1999-12-31T18:20:00.000000000', '1999-12-31T18:40:00.000000000',
           '1999-12-31T19:00:00.000000000', '1999-12-31T19:20:00.000000000',
           '1999-12-31T19:40:00.000000000', '1999-12-31T20:00:00.000000000',
           '1999-12-31T20:20:00.000000000', '1999-12-31T20:40:00.000000000',
           '1999-12-31T21:00:00.000000000', '1999-12-31T21:20:00.000000000',
           '1999-12-31T21:40:00.000000000', '1999-12-31T22:00:00.000000000',
           '1999-12-31T22:20:00.000000000', '1999-12-31T22:40:00.000000000',
           '1999-12-31T23:00:00.000000000', '1999-12-31T23:20:00.000000000',
           '1999-12-31T23:40:00.000000000', '2000-01-01T00:00:00.000000000',
           '2000-01-01T00:20:00.000000000', '2000-01-01T00:40:00.000000000',
           '2000-01-01T01:00:00.000000000', '2000-01-01T01:20:00.000000000',
           '2000-01-01T01:40:00.000000000', '2000-01-01T02:00:00.000000000',
           '2000-01-01T02:20:00.000000000', '2000-01-01T02:40:00.000000000',
           '2000-01-01T03:00:00.000000000', '2000-01-01T03:20:00.000000000',
           '2000-01-01T03:40:00.000000000', '2000-01-01T04:00:00.000000000',
           '2000-01-01T04:20:00.000000000', '2000-01-01T04:40:00.000000000',
           '2000-01-01T05:00:00.000000000', '2000-01-01T05:20:00.000000000',
           '2000-01-01T05:40:00.000000000', '2000-01-01T06:00:00.000000000',
           '2000-01-01T06:20:00.000000000', '2000-01-01T06:40:00.000000000',
           '2000-01-01T07:00:00.000000000', '2000-01-01T07:20:00.000000000',
           '2000-01-01T07:40:00.000000000', '2000-01-01T08:00:00.000000000',
           '2000-01-01T08:20:00.000000000', '2000-01-01T08:40:00.000000000',
           '2000-01-01T09:00:00.000000000', '2000-01-01T09:20:00.000000000',
           '2000-01-01T09:40:00.000000000', '2000-01-01T10:00:00.000000000',
           '2000-01-01T10:20:00.000000000', '2000-01-01T10:40:00.000000000',
           '2000-01-01T11:00:00.000000000', '2000-01-01T11:20:00.000000000',
           '2000-01-01T11:40:00.000000000', '2000-01-01T12:00:00.000000000',
           '2000-01-01T12:20:00.000000000', '2000-01-01T12:40:00.000000000',
           '2000-01-01T13:00:00.000000000', '2000-01-01T13:20:00.000000000',
           '2000-01-01T13:40:00.000000000', '2000-01-01T14:00:00.000000000',
           '2000-01-01T14:20:00.000000000', '2000-01-01T14:40:00.000000000',
           '2000-01-01T15:00:00.000000000', '2000-01-01T15:20:00.000000000',
           '2000-01-01T15:40:00.000000000', '2000-01-01T16:00:00.000000000',
           '2000-01-01T16:20:00.000000000', '2000-01-01T16:40:00.000000000',
           '2000-01-01T17:00:00.000000000', '2000-01-01T17:20:00.000000000',
           '2000-01-01T17:40:00.000000000', '2000-01-01T18:00:00.000000000',
           '2000-01-01T18:20:00.000000000', '2000-01-01T18:40:00.000000000',
           '2000-01-01T19:00:00.000000000', '2000-01-01T19:20:00.000000000',
           '2000-01-01T19:40:00.000000000', '2000-01-01T20:00:00.000000000',
           '2000-01-01T20:20:00.000000000', '2000-01-01T20:40:00.000000000',
           '2000-01-01T21:00:00.000000000', '2000-01-01T21:20:00.000000000',
           '2000-01-01T21:40:00.000000000', '2000-01-01T22:00:00.000000000',
           '2000-01-01T22:20:00.000000000', '2000-01-01T22:40:00.000000000',
           '2000-01-01T23:00:00.000000000', '2000-01-01T23:20:00.000000000',
           '2000-01-01T23:40:00.000000000', '2000-01-02T00:00:00.000000000',
           '2000-01-02T00:20:00.000000000', '2000-01-02T00:40:00.000000000',
           '2000-01-02T01:00:00.000000000', '2000-01-02T01:20:00.000000000',
           '2000-01-02T01:40:00.000000000', '2000-01-02T02:00:00.000000000',
           '2000-01-02T02:20:00.000000000', '2000-01-02T02:40:00.000000000',
           '2000-01-02T03:00:00.000000000', '2000-01-02T03:20:00.000000000',
           '2000-01-02T03:40:00.000000000', '2000-01-02T04:00:00.000000000',
           '2000-01-02T04:20:00.000000000', '2000-01-02T04:40:00.000000000',
           '2000-01-02T05:00:00.000000000', '2000-01-02T05:20:00.000000000',
           '2000-01-02T05:40:00.000000000', '2000-01-02T06:00:00.000000000',
           '2000-01-02T06:20:00.000000000', '2000-01-02T06:40:00.000000000',
           '2000-01-02T07:00:00.000000000', '2000-01-02T07:20:00.000000000',
           '2000-01-02T07:40:00.000000000', '2000-01-02T08:00:00.000000000',
           '2000-01-02T08:20:00.000000000', '2000-01-02T08:40:00.000000000',
           '2000-01-02T09:00:00.000000000', '2000-01-02T09:20:00.000000000',
           '2000-01-02T09:40:00.000000000', '2000-01-02T10:00:00.000000000',
           '2000-01-02T10:20:00.000000000', '2000-01-02T10:40:00.000000000',
           '2000-01-02T11:00:00.000000000', '2000-01-02T11:20:00.000000000',
           '2000-01-02T11:40:00.000000000', '2000-01-02T12:00:00.000000000',
           '2000-01-02T12:20:00.000000000', '2000-01-02T12:40:00.000000000',
           '2000-01-02T13:00:00.000000000', '2000-01-02T13:20:00.000000000',
           '2000-01-02T13:40:00.000000000', '2000-01-02T14:00:00.000000000',
           '2000-01-02T14:20:00.000000000', '2000-01-02T14:40:00.000000000',
           '2000-01-02T15:00:00.000000000', '2000-01-02T15:20:00.000000000',
           '2000-01-02T15:40:00.000000000', '2000-01-02T16:00:00.000000000'],
          dtype='datetime64[ns]')

  • SCHISM_hgrid_node_x
    (nSCHISM_hgrid_node)
    float64
    dask.array<chunksize=(2639,), meta=np.ndarray>

    axis :

    X

    location :

    node

    mesh :

    SCHISM_hgrid

    units :

    m

    standard_name :

    projection_x_coordinate

    Array Chunk Bytes 20.62 kiB 20.62 kiB Shape (2639,) (2639,) Dask graph 1 chunks in 2 graph layers Data type float64 numpy.ndarray

  • SCHISM_hgrid_node_y
    (nSCHISM_hgrid_node)
    float64
    dask.array<chunksize=(2639,), meta=np.ndarray>

    axis :

    Y

    location :

    node

    mesh :

    SCHISM_hgrid

    units :

    m

    standard_name :

    projection_y_coordinate

    Array Chunk Bytes 20.62 kiB 20.62 kiB Shape (2639,) (2639,) Dask graph 1 chunks in 2 graph layers Data type float64 numpy.ndarray

  • SCHISM_hgrid_face_x
    (nSCHISM_hgrid_face)
    float64
    dask.array<chunksize=(4636,), meta=np.ndarray>

    axis :

    X

    location :

    face

    mesh :

    SCHISM_hgrid

    units :

    m

    standard_name :

    projection_x_coordinate

    Array Chunk Bytes 36.22 kiB 36.22 kiB Shape (4636,) (4636,) Dask graph 1 chunks in 2 graph layers Data type float64 numpy.ndarray

  • SCHISM_hgrid_face_y
    (nSCHISM_hgrid_face)
    float64
    dask.array<chunksize=(4636,), meta=np.ndarray>

    axis :

    Y

    location :

    face

    mesh :

    SCHISM_hgrid

    units :

    m

    standard_name :

    projection_y_coordinate

    Array Chunk Bytes 36.22 kiB 36.22 kiB Shape (4636,) (4636,) Dask graph 1 chunks in 2 graph layers Data type float64 numpy.ndarray

  • SCHISM_hgrid_edge_x
    (nSCHISM_hgrid_edge)
    float64
    dask.array<chunksize=(7274,), meta=np.ndarray>

    axis :

    X

    location :

    edge

    mesh :

    SCHISM_hgrid

    units :

    m

    standard_name :

    projection_x_coordinate

    Array Chunk Bytes 56.83 kiB 56.83 kiB Shape (7274,) (7274,) Dask graph 1 chunks in 2 graph layers Data type float64 numpy.ndarray

  • SCHISM_hgrid_edge_y
    (nSCHISM_hgrid_edge)
    float64
    dask.array<chunksize=(7274,), meta=np.ndarray>

    axis :

    Y

    location :

    edge

    mesh :

    SCHISM_hgrid

    units :

    m

    standard_name :

    projection_y_coordinate

    Array Chunk Bytes 56.83 kiB 56.83 kiB Shape (7274,) (7274,) Dask graph 1 chunks in 2 graph layers Data type float64 numpy.ndarray

• Data variables: (12)
  • minimum_depth
    (one)
    float64
    dask.array<chunksize=(1,), meta=np.ndarray>

    units :

    m

    Array Chunk Bytes 8 B 8 B Shape (1,) (1,) Dask graph 1 chunks in 2 graph layers Data type float64 numpy.ndarray

  • SCHISM_hgrid
    (one)
    |S1
    dask.array<chunksize=(1,), meta=np.ndarray>

    long_name :

    Topology data of 2d unstructured mesh

    topology_dimension :

    2

    cf_role :

    mesh_topology

    node_coordinates :

    SCHISM_hgrid_node_x SCHISM_hgrid_node_y

    edge_coordinates :

    SCHISM_hgrid_edge_x SCHISM_hgrid_edge_y

    face_coordinates :

    SCHISM_hgrid_face_x SCHISM_hgrid_face_y

    edge_node_connectivity :

    SCHISM_hgrid_edge_nodes

    face_node_connectivity :

    SCHISM_hgrid_face_nodes

    Array Chunk Bytes 1 B 1 B Shape (1,) (1,) Dask graph 1 chunks in 2 graph layers Data type |S1 numpy.ndarray

  • depth
    (nSCHISM_hgrid_node)
    float32
    dask.array<chunksize=(2639,), meta=np.ndarray>

    units :

    m

    axis :

    Z

    positive :

    down

    location :

    node

    grid_mapping :

    crs

    mesh :

    SCHISM_hgrid

    Array Chunk Bytes 10.31 kiB 10.31 kiB Shape (2639,) (2639,) Dask graph 1 chunks in 2 graph layers Data type float32 numpy.ndarray

  • bottom_index_node
    (nSCHISM_hgrid_node)
    int32
    dask.array<chunksize=(2639,), meta=np.ndarray>

    location :

    node

    grid_mapping :

    crs

    mesh :

    SCHISM_hgrid

    Array Chunk Bytes 10.31 kiB 10.31 kiB Shape (2639,) (2639,) Dask graph 1 chunks in 2 graph layers Data type int32 numpy.ndarray

  • SCHISM_hgrid_face_nodes
    (nSCHISM_hgrid_face, nMaxSCHISM_hgrid_face_nodes)
    float64
    dask.array<chunksize=(4636, 4), meta=np.ndarray>

    start_index :

    1

    cf_role :

    face_node_connectivity

    Array Chunk Bytes 144.88 kiB 144.88 kiB Shape (4636, 4) (4636, 4) Dask graph 1 chunks in 2 graph layers Data type float64 numpy.ndarray

  • SCHISM_hgrid_edge_nodes
    (nSCHISM_hgrid_edge, two)
    float64
    dask.array<chunksize=(7274, 2), meta=np.ndarray>

    start_index :

    1

    cf_role :

    edge_node_connectivity

    Array Chunk Bytes 113.66 kiB 113.66 kiB Shape (7274, 2) (7274, 2) Dask graph 1 chunks in 2 graph layers Data type float64 numpy.ndarray

  • dryFlagNode
    (time, nSCHISM_hgrid_node)
    float32
    dask.array<chunksize=(72, 2639), meta=np.ndarray>

    i23d :

    1

    location :

    node

    grid_mapping :

    crs

    mesh :

    SCHISM_hgrid

    Array Chunk Bytes 1.45 MiB 742.22 kiB Shape (144, 2639) (72, 2639) Dask graph 2 chunks in 5 graph layers Data type float32 numpy.ndarray

  • elevation
    (time, nSCHISM_hgrid_node)
    float32
    dask.array<chunksize=(72, 2639), meta=np.ndarray>

    i23d :

    1

    location :

    node

    grid_mapping :

    crs

    mesh :

    SCHISM_hgrid

    Array Chunk Bytes 1.45 MiB 742.22 kiB Shape (144, 2639) (72, 2639) Dask graph 2 chunks in 5 graph layers Data type float32 numpy.ndarray

  • depthAverageVelX
    (time, nSCHISM_hgrid_node)
    float32
    dask.array<chunksize=(72, 2639), meta=np.ndarray>

    i23d :

    1

    location :

    node

    grid_mapping :

    crs

    mesh :

    SCHISM_hgrid

    Array Chunk Bytes 1.45 MiB 742.22 kiB Shape (144, 2639) (72, 2639) Dask graph 2 chunks in 5 graph layers Data type float32 numpy.ndarray

  • depthAverageVelY
    (time, nSCHISM_hgrid_node)
    float32
    dask.array<chunksize=(72, 2639), meta=np.ndarray>

    i23d :

    1

    location :

    node

    grid_mapping :

    crs

    mesh :

    SCHISM_hgrid

    Array Chunk Bytes 1.45 MiB 742.22 kiB Shape (144, 2639) (72, 2639) Dask graph 2 chunks in 5 graph layers Data type float32 numpy.ndarray

  • dryFlagElement
    (time, nSCHISM_hgrid_face)
    float32
    dask.array<chunksize=(72, 4636), meta=np.ndarray>

    i23d :

    4

    location :

    face

    grid_mapping :

    crs

    mesh :

    SCHISM_hgrid

    Array Chunk Bytes 2.55 MiB 1.27 MiB Shape (144, 4636) (72, 4636) Dask graph 2 chunks in 5 graph layers Data type float32 numpy.ndarray

  • dryFlagSide
    (time, nSCHISM_hgrid_edge)
    float32
    dask.array<chunksize=(72, 7274), meta=np.ndarray>

    i23d :

    7

    location :

    edge

    grid_mapping :

    crs

    mesh :

    SCHISM_hgrid

    Array Chunk Bytes 4.00 MiB 2.00 MiB Shape (144, 7274) (72, 7274) Dask graph 2 chunks in 5 graph layers Data type float32 numpy.ndarray

• Indexes: (1)
  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['1999-12-31 16:20:00', '1999-12-31 16:40:00',
                   '1999-12-31 17:00:00', '1999-12-31 17:20:00',
                   '1999-12-31 17:40:00', '1999-12-31 18:00:00',
                   '1999-12-31 18:20:00', '1999-12-31 18:40:00',
                   '1999-12-31 19:00:00', '1999-12-31 19:20:00',
                   ...
                   '2000-01-02 13:00:00', '2000-01-02 13:20:00',
                   '2000-01-02 13:40:00', '2000-01-02 14:00:00',
                   '2000-01-02 14:20:00', '2000-01-02 14:40:00',
                   '2000-01-02 15:00:00', '2000-01-02 15:20:00',
                   '2000-01-02 15:40:00', '2000-01-02 16:00:00'],
                  dtype='datetime64[ns]', name='time', length=144, freq=None))

• Attributes: (0)

Read mesh data of elements and nodes

This might not work for quad elements. Will have to deal with those as two non-overlapping triangles

smesh = schism_mesh.read_mesh('../tests/data/m1_hello_schism/hgrid.gr3')

dfelems = pd.DataFrame(smesh.elems,columns=[0,1,2])
#dfelems

dfnodes = pd.DataFrame(smesh.nodes, columns=['x','y','z'])
dfnodes.z = -dfnodes.z
#dfnodes

TriSurface plots from Plotly

Needed to add simplices from existing elements information of mesh rather than Delaunay triangulation of nodes

Adapted from https://anaconda.org/philippjfr/brain/notebook?version=2017.05.04.1924

The code below allows for the simplices already defined by elems to be used instead of doing a Delaunay triangulation (used from scipy as a way to calculate the simplices)

hv.extension('plotly')

import param

class TriSurface(hv.TriSurface):
    
    simplices = param.Array()

class TriSurfacePlot(hv.plotting.plotly.TriSurfacePlot):

    style_opts = ['cmap', 'plot_edges']

    def get_data(self, element, ranges, style, **kwargs):
        if element.simplices is None:
            return super(TriSurfacePlot, self).get_data(element, ranges, style, **kwargs)
        x, y, z = (element.dimension_values(i) for i in range(3))
        simplices = element.simplices
        return [dict(x=x, y=y, z=z, simplices=simplices)]
    
hv.Store.register({TriSurface: TriSurfacePlot}, 'plotly')

dfelev = dfnodes.copy()

#ds.elevation.min(), ds.elevation.max()

mesh_surface = TriSurface(dfelev, simplices = dfelems.values).opts(plot_edges=False, cmap='gray')
mesh_surface

Show mesh and water surface

The water surface simplices are derived from the nodes with Delaunay triangulation. Can we get the simplices from the information in the schism output files?

dfsurface = dfnodes.copy()

def show_surface(time=0):
    dfsurface.z = ds.elevation.values[time,:]
    water_surface = TriSurface(dfsurface, simplices = dfelems.values).opts(width=800, zlim=(-10,2))
    return water_surface.opts(plot_edges=False, cmap='kbc', clim=(0,2), colorbar=True)

show_surface(0)

Now combine both to animate both bathymetry and water surface

def show_combined(time):
    return show_surface(time)*mesh_surface

import panel as pn
pn.extension()

time_slider = pn.widgets.IntSlider(name='Time Index', start=0, end=len(ds.time)-1)

pn.Row(time_slider, 
       hv.DynamicMap(show_combined, streams={'time':time_slider}).opts(
           width=800, height=800)).servable(title='Water Surface Elevation Animation in 3D')#.show()